program Main_elast

    use omp_lib
    use Globals
    use ModuleElasticities
    use ModuleSolve_GLO_LOC
    use ModuleSolveEqm
    
    implicit none
    
    ! Local variable declarations
    real(8) :: param_in_vec(3)  ! vector of tax function parameters: input for function solving for equilibrium given taxes
    real(8) :: welfare          ! welfare associated with equilibrium
    real(8) :: ParamCalib(5), eqm_results(5)    ! calibrated parameters    

    ! Set tax function parameters
    open(11,  file = 'Input/eqm_results.txt',      status = 'unknown')
    read(11, *) eqm_results
    close(11)
    norm_version = 1            ! indicator for choice of normalization: 1) mean income; 2) median income; 3) mean income of the third quintile
    param_in_vec(1) = eqm_results(3)     ! gamma
    param_in_vec(2) = eqm_results(2)     ! lambda
    param_in_vec(3) = eqm_results(5)     ! rt
    omegat = 0d0                         ! omegat
    
    ! equilibrium guesses
    mt = eqm_results(4)         ! level parameter transfer function
    r = eqm_results(1)          ! interest rate

    ! Calibrated parameters
    open(11,  file = 'Input/ParamCalib.txt',      status = 'unknown')
    read(11, *) ParamCalib
    close(11)
    B           = ParamCalib(2)             ! labor disutility parameter
    beta        = ParamCalib(1)             ! discount factor
    G           = ParamCalib(3)             ! government spending
    Debt        = ParamCalib(4)             ! government debt
    ymean       = ParamCalib(5)             ! mean income
    ymeanCalib  = ParamCalib(5)             ! mean income
    
    ! Solve for equilibrium given tax function parameters
    epsilon = 0d0                           ! no change to tax system
    index_epsilon = 1                       ! initial steady state
    welfare = welfare_eqm(param_in_vec,1)   ! solve

    ! First call of elasticities to store some variables
    call elasticities

    ! Solve for policies given marginal change to taxes
    epsilon = 0.01d0                        ! small change to tax system
    index_epsilon = 2                       ! changed tax system
    call Solve_GLO_LOC                      ! solve HH problem

    ! Second call of elasticities to compute elasticities
    call elasticities
    
end program Main_elast
